/*
* INSANE - Interactive Structural Analysis Environment
*
* Copyright (C) 2003-2005
* Universidade Federal de Minas Gerais
* Escola de Engenharia
* Departamento de Engenharia de Estruturas
* 
* Author's email :    insane@dees.ufmg.br
* Author's website :  http://www.dees.ufmg.br/insane
* 
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License
* as published by the Free Software Foundation; either version 2
* of the License, or any later version.
* 
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
* 
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
*/
package br.ufmg.dees.insane.util;

import java.io.Serializable;
import java.text.DecimalFormat;
import java.util.ArrayList;

/**A double precision, general, dynamically-resizable, two-dimensional matrix class.
* Row and column numbering begins with zero.
* The representation is row major.
* @author Fonseca, Flavio & Pitangueira, Roque
* @since May 2004
*/
public class IMatrix implements Serializable, Cloneable {
	
	/**
     * Comment for <code>serialVersionUID</code>
     */
    private static final long serialVersionUID = 1L;
    int numRow;
	int numCol;
	double[][] matrix;
 //*****************************************************************************
	
	/**Constructs a new IMatrix and copies the initial values from the parameter matrix.
	*@param matrix the source of the initial values of the new IMatrix
	*/
	public IMatrix(IMatrix matrix) {
		numRow = matrix.getNumRow();
		numCol = matrix.getNumCol();
		this.matrix = new double[numRow][numCol];
		for (int i=0; i<numRow; i++) {
			for (int j=0; j<numCol; j++) {
				this.matrix[i][j] = matrix.getElement(i, j);
			}
		}
	}
    
    /**Constructs a new IMatrix initialized to the values in the matrix array.
    *@param matrix a 2D array that specifies a matrix
    */
    public IMatrix(double[][] matrix) {
        numRow = matrix.length;
        numCol = matrix.length;
        // So para matriz quadrada
        this.matrix = new double[numRow][numCol];
        for (int i=0; i<numRow; i++) {
            for (int j=0; j<numCol; j++) {
                this.matrix[i][j] = matrix[i][j];
            }
        }
    }
	
//*****************************************************************************
	
	/**Constructs an nRow by nCol identity matrix.
	*@param nRow number of rows in this matrix.
	*@param nCol number of columns in this matrix.
	*/
	public IMatrix(int nRow, int nCol) {
		numRow = nRow;
		numCol = nCol;
		matrix = new double[numRow][numCol];
		this.setIdentity();
	}
	
//*****************************************************************************
	
	/**Constructs an nRow by nCol matrix initialized to the values in the matrix array.
	*@param nRow number of rows in this matrix.
	*@param nCol number of columns in this matrix.
	*@param matrix a 1D array that specifies a matrix in row major fashion
	*/
	public IMatrix(int nRow, int nCol, double[] matrix) {
		numRow = nRow;
		numCol = nCol;
		this.matrix = new double[numRow][numCol];
		int k = 0;
		for (int i=0; i<numRow; i++) {
			for (int j=0; j<numCol; j++) {
				this.matrix[i][j] = matrix[k];
				k++;
			}
		}
		this.setIdentity();
	}
	
//*****************************************************************************
	
	/**Sets the value of this matrix to sum of itself and matrix m1.
	*@param m1 the other matrix
	*/
	public final void add(IMatrix m1) throws MismatchedSizeException{
		
		if (m1.getNumRow()!=this.numRow || m1.getNumCol()!=this.numCol) {
			throw (new MismatchedSizeException());
		}
		
		for (int i=0; i<this.numRow; i++) {
			for (int j=0; j<this.numCol; j++) {
				this.matrix[i][j] = this.matrix[i][j] + m1.getElement(i, j);
			}
		}
	}
	
//*****************************************************************************
	
	/**Sets the value of this matrix to the matrix sum of matrices m1 and m2.
	*@param m1 the first matrix
	*@param m2 the second matrix
	*/
	public final void add(IMatrix m1, IMatrix m2) throws MismatchedSizeException{
		
		if (m1.getNumRow()!=this.numRow || m1.getNumCol()!=this.numCol) {
			throw (new MismatchedSizeException());
		}
		if (m1.getNumRow()!=m2.getNumRow() || m1.getNumCol()!=m2.getNumCol()) {
			throw (new MismatchedSizeException());
		}
		
		for (int i=0; i<this.numRow; i++) {
			for (int j=0; j<this.numCol; j++) {
				this.matrix[i][j] = m1.getElement(i, j) + m2.getElement(i, j);
			}
		}
	}
	
//*****************************************************************************
	
	/**Creates a new object of the same class as this object.*/
	public Object clone() {
		IMatrix aux = new IMatrix(this);
		return aux;
	}
	
//*****************************************************************************
	
	/**Returns true if all of the data members of IMatrix m1 are equal to the corresponding data members in this IMatrix.
	*@param m1 The matrix with which the comparison is made
	*@return true or false
	*/
	public boolean equals(IMatrix m1) {
		if (this.numRow != m1.getNumRow())
			return false;
		if (this.numCol != m1.getNumCol())
			return false;
		for (int i=0; i<this.numRow; i++)
			for (int j=0; j<this.numCol; j++)
				if (this.matrix[i][j] != m1.getElement(i,j))
					return false;
		return true;
	}
	
//*****************************************************************************
	
	/**Returns true if the Object o1 is of type IMatrix and all of the data members of o1 are equal to the corresponding data members in this IMatrix.
	*@param o1 The object with which the comparison is made
	*@return true or false
	*/
	public boolean equals(Object o1) {
		if (o1 instanceof IMatrix) {
			if (this.numRow != ((IMatrix)o1).getNumRow())
				return false;
			if (this.numCol != ((IMatrix)o1).getNumCol())
				return false;
				for (int i=0; i<this.numRow; i++)
					for (int j=0; j<this.numCol; j++)
						if (this.matrix[i][j] != ((IMatrix)o1).getElement(i,j))
							return false;
		}
		return true;
	}
	
//*****************************************************************************
	
	/**Places the values of the specified column into the array parameter.
	*@param col the target column number
	*@param	array the array into which the column values will be placed
	*/
	public final void getColumn(int col, double[] array) {
		for (int i=0; i<numRow; i++) {
			array[i] = this.matrix[i][col];
		}
	}
	
//*****************************************************************************
	
	/**Places the values of the specified column into the vector parameter.
	*@param col the target column number
	*@param vector the vector into which the column values will be placed
	*/
	public final void getColumn(int col, IVector vector) {
		for (int i=0; i<numRow; i++) {
			vector.setElement(i, this.matrix[i][col]);
		}
	}
	
//*****************************************************************************
	
	/**Places the values of the specified row into the array parameter.
	*@param row the target column number
	*@param	array the array into which the row values will be placed
	*/
	public final void getRow(int row, double[] array) {
		for (int j=0; j<numCol; j++) {
			array[j] = this.matrix[row][j];
		}
	}
	
//*****************************************************************************
	
	/**Places the values of the specified row into the vector parameter.
	*@param row the target column number
	*@param vector the vector into which the row values will be placed
	*/
	public final void getRow(int row, IVector vector) {
		for (int j=0; j<numCol; j++) {
			vector.setElement(j, this.matrix[row][j]);
		}
	}
	
//*****************************************************************************
	
	/**Places the values of the specified row into the vector parameter.
	*@param row the target column number
	*@param vector the vector into which the row values will be placed
	*/
	public final ArrayList<Double> getRow(int row) {
		ArrayList<Double> a = new ArrayList<Double>(this.getNumCol());
		for (int j=0; j<numCol; j++) {
			a.add(j, this.matrix[row][j]);
		}
		return a;
	}
	
//*****************************************************************************
	
	/**Retrieves the value at the specified row and column of this matrix.
	*@param row the row number to be retrieved (zero indexed)
	*@param column the column number to be retrieved (zero indexed)
	*@return the value at the indexed element
	*/
	public final double getElement(int row, int column) {
		return matrix[row][column];
	}
	
//*****************************************************************************
	
	/**Returns the number of colmuns in this matrix.
	*@return number of columns in this matrix
	*/
	public final int getNumCol() {
		return numCol;
	}
	
//*****************************************************************************
	
	/**Returns the number of rows in this matrix.
	*@return number of rows in this matrix
	*/
	public final int getNumRow() {
		return numRow;
	}
	
//*****************************************************************************
	
	/**Subtracts this matrix from the identity matrix and puts the values back into this (this = I - this).*/
	public final void identityMinus() {
		IMatrix aux = new IMatrix(this.numRow, this.numCol);
		aux.setIdentity();
		this.sub(aux, this);
	}
	
//*****************************************************************************
/*
void invert() 
Inverts this matrix in place.
*/
//*****************************************************************************
/*
void invert(IMatrix m1) 
Inverts matrix m1 and places the new values into this matrix.
*/
//*****************************************************************************
	
	/** Calculates the determinant of this matrix using LU decomposition.
	*@return the determinant of this matrix.
	*/
	public double determinantLU() throws MismatchedSizeException{
		
		if (this.numRow!=this.numCol){
			throw (new MismatchedSizeException("IMatrix.determinant() demands a square matrix."));
		}
		
		IMatrix LU = new IMatrix(this);
		IVector permutation = new IVector(this.numCol);
		
		for (int i=0; i<this.numRow; i++)
			for (int j=0; j<this.numCol; j++)
				LU.setElement(i, j, this.matrix[i][j]);
		
		double det;
		int n = LU.getNumRow()-1;
		
		for (int i=0; i<=n; i++)
			permutation.setElement(i, i);
		det = 1;
		
		for (int j=0; j<=(n-1); j++) {
			int p = 0;
			p = j;
			double amax = Math.abs(LU.getElement(j,j));
			for (int k=j+1; k<=n; k++) {
				if (Math.abs(LU.getElement(k,j)) > amax) {
					amax = Math.abs(LU.getElement(k,j));
					p = k;
				}
			}
			if (p!=j) {
				double t = 0;
				for (int k=0; k<=n; k++) {
					t = LU.getElement(j, k);
					LU.setElement(j,k, LU.getElement(p, k));
					LU.setElement(p,k,t);
				}
				t = permutation.getElement(j);
				permutation.setElement(j, permutation.getElement(p));
				permutation.setElement(p, t);
				det = -det;
			}
			det = det * LU.getElement(j,j);
			if (Math.abs(LU.getElement(j,j)) != 0) {
				double r = 1 / LU.getElement(j,j);
				for (int i=j+1; i<=n; i++) {
					double m = LU.getElement(i,j) * r;
					LU.setElement(i,j,m);
					for (int k=j+1; k<=n; k++) {
						LU.setElement(i, k, LU.getElement(i,k)-m*LU.getElement(j,k));
					}
				}
			} else {
				det = 0;
			}
		}
		det = det * LU.getElement(n,n);
		return det;
	}
	
//*****************************************************************************
	
	/**LU Decomposition: this matrix must be a square matrix and the LU IMatrix parameter must be the same size as this matrix.
	* The matrix LU will be overwritten as the combination of a lower diagonal and upper diagonal matrix decompostion of this matrix;
	* the diagonal elements of L (unity) are not stored.
	* The IVector parameter records the row permutation effected by the partial pivoting, and is used as a parameter to the Iector method LUDBackSolve to solve sets of linear equations.
	*@param LU The matrix into which the lower and upper decompositions will be placed.
	*@param permutation The row permutation effected by the partial pivoting 
	*/
	public final void LUD(IMatrix LU, IVector permutation) throws MismatchedSizeException{
		
		if (this.numRow!=this.numCol){
			throw (new MismatchedSizeException("IMatrix.LUD() demands a square matrix."));
		}
		if (this.numRow!=LU.getNumRow() || this.numCol!=LU.getNumCol() || this.numRow!=permutation.getSize()) {
			throw (new MismatchedSizeException());
		}
		
		for (int i=0; i<this.numRow; i++)
			for (int j=0; j<this.numCol; j++)
				LU.setElement(i, j, this.matrix[i][j]);
		
		double det;
		int n = LU.getNumRow()-1;
		
		for (int i=0; i<=n; i++)
			permutation.setElement(i, i);
		det = 1;
		
		for (int j=0; j<=(n-1); j++) {
			int p = 0;
			p = j;
			double amax = Math.abs(LU.getElement(j,j));
			for (int k=j+1; k<=n; k++) {
				if (Math.abs(LU.getElement(k,j)) > amax) {
					amax = Math.abs(LU.getElement(k,j));
					p = k;
				}
			}
			if (p!=j) {
				double t = 0;
				for (int k=0; k<=n; k++) {
					t = LU.getElement(j, k);
					LU.setElement(j,k, LU.getElement(p, k));
					LU.setElement(p,k,t);
				}
				t = permutation.getElement(j);
				permutation.setElement(j, permutation.getElement(p));
				permutation.setElement(p, t);
				det = -det;
			}
			det = det * LU.getElement(j,j);
			if (Math.abs(LU.getElement(j,j)) != 0) {
				double r = 1 / LU.getElement(j,j);
				for (int i=j+1; i<=n; i++) {
					double m = LU.getElement(i,j) * r;
					LU.setElement(i,j,m);
					for (int k=j+1; k<=n; k++) {
						LU.setElement(i, k, LU.getElement(i,k)-m*LU.getElement(j,k));
					}
				}
			} else {
				det = 0;
			}
		}
		det = det * LU.getElement(n,n);
	}
	
//*****************************************************************************
	
	/**Cholesky Decomposition: this matrix must be a symetric square matrix and the L IMatrix parameter must be the same size as this matrix.
	* The matrix L will be overwritten as the lower diagonal matrix obtained from the decompostion of this matrix;
	* The upper diagonal matrix is the transpose of the lower diagonal one.
	*@param L The matrix into which the lower decomposition will be placed.
	*/
	public final void choleskyD(IMatrix L) throws MismatchedSizeException, SingularMatrixException{
		
		if (this.numRow!=this.numCol){
			throw (new MismatchedSizeException("IMatrix.choleskyD() demands a square matrix."));
		}
		if (this.numRow!=L.getNumRow() || this.numCol!=L.getNumCol()) {
			throw (new MismatchedSizeException());
		}
		
		L.setZero();
		
		double det = 1;
		int n = L.getNumRow()-1;
		
		for (int j=0; j<=n; j++) {
			double sum = 0;
			for (int k=0; k<=(j-1); k++)
				sum = sum + Math.pow(L.getElement(j,k),2);
			double t = this.getElement(j,j) - sum;
			det = det * t;
			boolean error = false;
			double r;
			if (t<=0)
				error = true;
			if (error) {
				throw new SingularMatrixException();
			}
			else {
				L.setElement(j,j, Math.sqrt(t));
				r = 1 / L.getElement(j,j);
			}
			for (int i=j+1; i<=n; i++) {
				sum = 0;
				for (int k=0; k<=j-1; k++)
					sum = sum + L.getElement(i,k)*L.getElement(j,k);
				L.setElement(i,j, (this.getElement(i,j)-sum)*r);
			}
		}
		
		//Avoiding residual numbers problem
		double higher = L.getElement(0,0);
		double lower = L.getElement(0,0);
		double division = 0;
		for (int i=0; i<numRow; i++) {
			if (L.getElement(i,i) > higher)
				higher = L.getElement(i,i);
			if (L.getElement(i,i) < lower)
				lower = L.getElement(i,i);
		}
		division = lower / higher;
		if (division < 1e-5)
			throw new SingularMatrixException();
		
	}
	
//*****************************************************************************
	
	/** Calculates the determinant of this matrix using Cholesky decomposition.
	*@return the determinant of this matrix.
	*/
	public final double determinantCho() throws MismatchedSizeException, SingularMatrixException{
		
		if (this.numRow!=this.numCol){
			throw (new MismatchedSizeException("IMatrix.determinantCho() demands a square matrix."));
		}
		
		IMatrix L = new IMatrix(this.getNumRow(), this.getNumCol());
		
		double det = 1;
		int n = L.getNumRow()-1;
		
		for (int j=0; j<=n; j++) {
			double sum = 0;
			for (int k=0; k<=(j-1); k++)
				sum = sum + Math.pow(L.getElement(j,k),2);
			double t = this.getElement(j,j) - sum;
			det = det * t;
			boolean error = false;
			double r;
			if (t<=0)
				error = true;
			if (error)
				throw new SingularMatrixException();
			else {
				L.setElement(j,j, Math.sqrt(t));
				r = 1 / L.getElement(j,j);
			}
			for (int i=j+1; i<=n; i++) {
				sum = 0;
				for (int k=0; k<=j-1; k++)
					sum = sum + L.getElement(i,k)*L.getElement(j,k);
				L.setElement(i,j, (this.getElement(i,j)-sum)*r);
			}
		}
		
		return det;
	}
	
//*****************************************************************************
	
	/**Sets the value of this matrix to the result of multiplying itself with vector v1 (this = this * v1).
	*@param v1 the vector
	*/
	public final IVector mul(IVector v1) throws MismatchedSizeException{
		if (this.numCol!=v1.getSize()) {
			throw (new MismatchedSizeException());
		}
		IVector aux = new IVector(this.numRow);
		for (int i=0; i<this.numRow; i++) {
				for (int k=0; k<this.numCol; k++) {
						aux.setElement(i,  aux.getElement(i) + this.matrix[i][k] * v1.getElement(k));
				}
		}
		return aux;
	}
	
	/**Sets the value of this matrix to the result of multiplying itself with matrix m1 (this = this * m1).
	*@param m1 the other matrix
	*/
	public final void mul(IMatrix m1) throws MismatchedSizeException{
		if (this.numCol!=m1.getNumRow()) {
			throw (new MismatchedSizeException());
		}
		IMatrix aux = new IMatrix(this.numRow, m1.getNumCol());
		aux.setZero();
		for (int i=0; i<this.numRow; i++) {
		    for (int j=0; j<m1.getNumCol(); j++) {
		        for (int k=0; k<this.numCol; k++) {
		            aux.setElement(i, j,  aux.getElement(i,j) + this.matrix[i][k] * m1.getElement(k, j));
		        }
		    }
		}
		this.setSize(this.numRow, m1.getNumCol());
		for (int i=0; i<this.numRow; i++)
			for (int j=0; j<m1.getNumCol(); j++)
				this.matrix[i][j] = aux.getElement(i,j);
	}
	
//*****************************************************************************
	
	/**Sets the value of this matrix to the result of multiplying the two argument matrices together (this = m1 * m2).
	*@param m1 the first matrix
	*@param m2 the second matrix
	*/
	public final void mul(IMatrix m1, IMatrix m2) throws MismatchedSizeException{
		
		if (m1.getNumCol()!=m2.getNumRow() || m1.getNumRow()!=this.numRow || m2.getNumCol()!=this.numCol) {
			throw (new MismatchedSizeException());
		}
		
		this.setZero();
		for (int i=0; i<m1.getNumRow(); i++)
			for (int j=0; j<m2.getNumCol(); j++)
				for (int k=0; k<m1.getNumCol(); k++)
						this.matrix[i][j] += m1.getElement(i, k) * m2.getElement(k, j);
	}
	
//*****************************************************************************
	
	/**Computes the outer product of the two vectors;
	* multiplies the the first vector by the transpose of the second vector and places the matrix result into this matrix.
	*@param v1 the first vector, treated as a row vector
	*@param v2 the second vector, treated as a column vector
	*/
	public final void mul(IVector v1, IVector v2) throws MismatchedSizeException{
		if (v1.getSize()!=v2.getSize() || v1.getSize()!=this.numRow || v2.getSize()!=this.numCol) {
				throw (new MismatchedSizeException());
		}
		
		for (int i=0; i<this.numRow; i++)
			for (int j=0; j<this.numCol; j++)
				this.matrix[i][j] = v1.getElement(i)*v2.getElement(j);
		
	}
	
//*****************************************************************************
	
	/**Multiplies the transpose of matrix m1 times the transpose of matrix m2, and places the result into this.
	*@param m1 The matrix on the left hand side of the multiplication
	*@param m2 The matrix on the right hand side of the multiplication
	*/
	public final void mulTransposeBoth(IMatrix m1, IMatrix m2) throws MismatchedSizeException{
		if (m1.getNumRow()!=m2.getNumCol() || m1.getNumCol()!=this.numRow || m2.getNumRow()!=this.numCol) {
				throw (new MismatchedSizeException());
		}
		
		this.setZero();
		IMatrix m1a = new IMatrix(m1);
		IMatrix m2a = new IMatrix(m2);
		m1a.transpose();
		m2a.transpose();
		this.mul(m1a, m2a);
	}
	
//*****************************************************************************
	
	/**Multiplies the transpose of matrix m1 times matrix m2, and places the result into this.
	*@param m1 The matrix on the left hand side of the multiplication
	*@param m2 The matrix on the right hand side of the multiplication
	*/
	public final void mulTransposeLeft(IMatrix m1, IMatrix m2) throws MismatchedSizeException{
		
		if (m1.getNumRow()!=m2.getNumRow() || m1.getNumCol()!=this.numRow || m2.getNumCol()!=this.numCol) {
			throw (new MismatchedSizeException());
		}
		
		this.setZero();
		IMatrix m1a = new IMatrix(m1);
		m1a.transpose();
		this.mul(m1a, m2);
	}
	
//*****************************************************************************
	
	/**Multiplies matrix m1 times the transpose of matrix m2, and places the result into this.
	*@param m1 The matrix on the left hand side of the multiplication
	*@param m2 The matrix on the right hand side of the multiplication
	*/
	public final void mulTransposeRight(IMatrix m1, IMatrix m2) throws MismatchedSizeException{
		if (m1.getNumCol()!=m2.getNumCol() || m1.getNumRow()!=this.numRow || m2.getNumRow()!=this.numCol) {
			throw (new MismatchedSizeException());
		}
		
		this.setZero();
		IMatrix m2a = new IMatrix(m2);
		m2a.transpose();
		this.mul(m1, m2a);
	}
	
//*****************************************************************************
	
	/**Negates the value of this matrix: this = -this.*/
	public final void negate() {
		for (int i=0; i<this.numRow; i++) {
			for (int j=0; j<this.numCol; j++) {
				this.matrix[i][j] = this.matrix[i][j] * -1;
			}
		}
	}
	
//*****************************************************************************
	
	/**Sets the value of this matrix equal to the negation of the IMatrix parameter.
	*@param m1 The source matrix
	*/
	public final void negate(IMatrix m1) {
		if (m1.getNumRow()!=this.numRow || m1.getNumCol()!=this.numCol) {
			throw (new MismatchedSizeException());
		}
		for (int i=0; i<this.numRow; i++) {
			for (int j=0; j<this.numCol; j++) {
				this.matrix[i][j] = m1.getElement(i,j) * -1;
			}
		}
	}
	
//*****************************************************************************
	
	/**Sets the value of this matrix to the values found in the array parameter.
	* The values are copied in one row at a time, in row major fashion.
	* The array should be at least equal in length to the number of matrix rows times the number of matrix columns in this matrix.
	*@param matrix the row major source array
	*/
	public final void set(double[] matrix) throws MismatchedSizeException{
		if (matrix.length < (this.numRow*this.numCol)){
			throw (new MismatchedSizeException());
		}
		int k = 0;
		for (int i=0; i<this.numRow; i++) {
			for (int j=0; j<this.numCol; j++) {
				this.matrix[i][j] = matrix[k];
				k++;
			}
		}
	}
	
//*****************************************************************************
	
	/**Sets the value of this matrix to the values found in matrix m1.
	*@param m1 the source matrix
	*/
	public final void set(IMatrix m1) throws MismatchedSizeException{
		if (m1.getNumRow()!=this.numRow || m1.getNumCol()!=this.numCol) {
			throw (new MismatchedSizeException());
		}
		for (int i=0; i<this.numRow; i++) {
			for (int j=0; j<this.numCol; j++) {
				this.matrix[i][j] = m1.getElement(i,j);
			}
		}
	}
	
//*****************************************************************************
	
	/**Copy the values from the array into the specified column of this matrix.
	*@param col the column of this matrix into which the array values will be copied
	*@param array the source array
	*/
	public final void setColumn(int col, double[] array) {
		for (int i=0; i<numRow; i++) {
			this.matrix[i][col] = array[i];
		}
	}
	
//*****************************************************************************
	
	/**Copy the values from the vector into the specified column of this matrix.
	*@param col the column of this matrix into which the array values will be copied
	*@param vector the source vector
	*/
	public final void setColumn(int col, IVector vector) {
		for (int i=0; i<numRow; i++) {
			this.matrix[i][col] = vector.getElement(i);
		}
	}
	
//*****************************************************************************
	
	/**Modifies the value at the specified row and column of this matrix.
	*@param row the row number to be modified (zero indexed)
	*@param column the column number to be modified (zero indexed)
	*@param value the new matrix element value
	*/
	public final void setElement(int row, int column, double value) {
		this.matrix[row][column] = value;
	}
	
//*****************************************************************************
	
	/**Sets this IMatrix to the identity matrix.*/
	public final void setIdentity() {
		for (int i=0; i<this.numRow; i++) {
			for (int j=0; j<this.numCol; j++) {
				if (i==j)
					this.matrix[i][j]=1;
				else
					this.matrix[i][j] = 0;
			}
		}
	}
	
//*****************************************************************************
	
	/**Copy the values from the array into the specified row of this matrix.
	*@param row the row of this matrix into which the array values will be copied.
	*@param array the source array
	*/
	public final void setRow(int row, double[] array) {
		for (int j=0; j<numCol; j++) {
			this.matrix[row][j] = array[j];
		}
	}
	
//*****************************************************************************
	
	/**Copy the values from the vector into the specified row of this matrix.
	*@param row the row of this matrix into which the array values will be copied
	*@param vector the source vector
	*/
	public final void setRow(int row, IVector vector) {
		for (int j=0; j<numCol; j++) {
			this.matrix[row][j] = vector.getElement(j);
		}
	}
	
//*****************************************************************************
	
	/**Sets this matrix to a uniform scale matrix; all of the values are reset.
	*@param scale The new scale value
	*/
	public final void setScale(double scale) {
		for (int i=0; i<this.numRow; i++) {
			for (int j=0; j<this.numCol; j++) {
				this.matrix[i][j] = this.matrix[i][j] * scale;
			}
		}
	}
	
//*****************************************************************************
	
	/**Changes the size of this matrix dynamically.
	* If the size is increased no data values will be lost.
	* If the size is decreased, only those data values whose matrix positions were eliminated will be lost.
	*@param nRow number of desired rows in this matrix
	*@param nCol number of desired columns in this matrix
	*/
	public void setSize(int nRow, int nCol) {
	    double[][] aux = new double[nRow][nCol];
	    
	    if (nRow<this.numRow && nCol<this.numCol) {
			for (int i=0; i<nRow; i++) {
		        for (int j=0; j<nCol; j++) {
					        aux[i][j]=this.matrix[i][j];
		                }
			}
	    }
	    if (nRow>=this.numRow && nCol<this.numCol) {
			for (int i=0; i<this.numRow; i++) {
		        for (int j=0; j<nCol; j++) {
					        aux[i][j]=this.matrix[i][j];
		                }
			}
	    }
	    if (nRow<this.numRow && nCol>=this.numCol) {
			for (int i=0; i<nRow; i++) {
		        for (int j=0; j<this.numCol; j++) {
					        aux[i][j]=this.matrix[i][j];
		                }
			}
	    }
	    if (nRow>=this.numRow && nCol>=this.numCol) {
			for (int i=0; i<this.numRow; i++) {
		        for (int j=0; j<this.numCol; j++) {
					        aux[i][j]=this.matrix[i][j];
		                }
			}
	    }
	    this.matrix= new double[nRow][nCol];
	    this.numCol=nCol;
	    this.numRow=nRow;
	    this.matrix=aux;
	}
	
//*****************************************************************************
	
	/**Sets all the values in this matrix to zero.*/
	public final void setZero(){
		for (int i=0; i<this.numRow; i++) {
			for (int j=0; j<this.numCol; j++) {
				this.matrix[i][j] = 0;
			}
		}
	}
	
//*****************************************************************************
	
	/**Sets the value of this matrix to the matrix difference of itself and matrix m1 (this = this - m1).
	*@param m1 - the other matrix
	*/
	public final void sub(IMatrix m1) throws MismatchedSizeException{
		if (m1.getNumRow()!=this.numRow || m1.getNumCol()!=this.numCol) {
			throw (new MismatchedSizeException());
		}
		for (int i=0; i<this.numRow; i++) {
			for (int j=0; j<this.numCol; j++) {
				this.matrix[i][j] = this.matrix[i][j] - m1.getElement(i, j);
			}
		}
	}
	
//*****************************************************************************
	
	/**Sets the value of this matrix to the matrix difference of matrices m1 and m2 (this = m1 - m2).
	*@param m1 the first matrix
	*@param m2 the second matrix
	*/
	public final void sub(IMatrix m1, IMatrix m2) throws MismatchedSizeException{
		if (m1.getNumRow()!=this.numRow || m1.getNumCol()!=this.numCol) {
			throw (new MismatchedSizeException());
		}
		if (m1.getNumRow()!=m2.getNumRow() || m1.getNumCol()!=m2.getNumCol()) {
			throw (new MismatchedSizeException());
		}
		for (int i=0; i<this.numRow; i++) {
			for (int j=0; j<this.numCol; j++) {
				this.matrix[i][j] = m1.getElement(i, j) - m2.getElement(i, j);
			}
		}
	}
	
//*****************************************************************************
	
	/**Transposes this matrix in place.*/
	public final void transpose() {
		IMatrix aux = new IMatrix(this.numRow, this.numCol);
		for (int i=0; i<this.numRow; i++)
			for (int j=0; j<this.numCol; j++)
				aux.setElement(i, j, this.getElement(i,j));
		this.transpose(aux);
	}
	
//*****************************************************************************
	
	/**Places the matrix values of the transpose of matrix m1 into this matrix.
	*@param m1 the matrix to be transposed (but not modified)
	*/
	public final void transpose(IMatrix m1) {
		this.setSize(m1.getNumCol(), m1.getNumRow());
		for (int i=0; i<m1.getNumRow(); i++) {
			for (int j=0; j<m1.getNumCol(); j++) {
				this.matrix[j][i]=m1.getElement(i,j);
			}
		}
	}
	
//*****************************************************************************
	
	/**Copies a sub-matrix derived from this matrix into the target matrix. The upper left of the sub-matrix is located at (rowSource, colSource); the lower right of the sub-matrix is located at (lastRowSource,lastColSource). The sub-matrix is copied into the the target matrix starting at (rowDest, colDest).
	*@param rowSource the top-most row of the sub-matrix
	*@param colSource the left-most column of the sub-matrix
	*@param numRow the number of rows in the sub-matrix
	*@param numCol the number of columns in the sub-matrix
	*@param rowDest the top-most row of the position of the copied sub-matrix within the target matrix
	*@param colDest the left-most column of the position of the copied sub-matrix within the target matrix
	*@param target the matrix into which the sub-matrix will be copied
	*/
	public final void copySubMatrix(int rowSource,int colSource, int numRow, int numCol, int rowDest, int colDest, IMatrix target) {
		
		for (int i=0; i<numRow; i++)
			for (int j=0; j<numCol; j++)
				target.setElement(rowDest+i, colDest+j, this.getElement(rowSource+i, colSource+j));
		
	}
	
//*****************************************************************************
	
	/**Places the column matrix values into a vector.
	*@return The vector containing the column matrix values (but not modified)
	*/
	public final IVector toVector() {
		IVector v = new IVector(this.getNumRow());
		for (int i=0; i<v.getSize(); i++) {
		    v.setElement(i,this.getElement(i,0));
		}
		return v;
	}
	
	/**Places the matrix values into a array of doubles.
	*@return The array of doubles containing the matrix values (but not modified)
	*/
	public final double[][] toDoubleArray() {
		double [][] d = new double[this.getNumRow()][this.getNumCol()];
		for (int i=0; i<this.getNumRow(); i++)
		    for (int j=0; j<this.getNumCol(); j++)
		{
		    d[i][j] = this.getElement(i,j);
		}
		return d;
	}

    /**Scales this matrix by the scale factor s.
    *@param s the scalar value
    */
    public final void scale(double s) {
        for (int i=0; i<this.getNumRow(); i++)
            for (int j=0; j<this.getNumCol(); j++){
                this.matrix[i][j] = this.matrix[i][j] * s;
        }
    }
    
    /** Returns the inverse of any square matrix
     * @return The inverse matrix of this matrix
     */
    public final IMatrix inverse(){
      if (this.numRow!=this.numCol){
          throw (new MismatchedSizeException("IMatrix.inverse() demands a square matrix."));
      }
      IMatrix inverse=new IMatrix(this.numRow,this.numCol);      
      IMatrix matL=new IMatrix(this.numRow,this.numCol);
      IMatrix matU=new IMatrix(this.numRow,this.numCol);
      this.genericL(matL,matU,1.0);      
      IMatrix matU1 = matU.inverseTriangleUpper(); 
      IMatrix matL1 = matL.inverseTriangleLower(); 
      inverse.mul(matU1,matL1);
      return inverse;
    }
   
    /** Returns the inverse of this triangle lower matrix       
     * @return The inverse matrix 
     */
    public final IMatrix inverseTriangleLower(){
      if (this.numRow!=this.numCol){
          throw (new MismatchedSizeException("IMatrix.inverseTriangularInferior() demands a square matrix."));
      }
      int n=this.numRow;
      IMatrix inverse=new IMatrix(n,n);
      for (int i=0;i<n;i++){
        for (int j=0;j<n;j++){
          if (i==j) inverse.setElement(i,j,1/this.matrix[i][i]);
          else if (i<j) inverse.setElement(i,j,0);
          else if (i>j) {
            double soma=0;
            for (int k=j;k<i;k++) soma+=this.matrix[i][k]*inverse.getElement(k,j);
            inverse.setElement(i,j,-soma/this.matrix[i][i]);
          }
        }
      }
      return inverse;
    }
    
    /** Returns the inverse of this triangle upper matrix       
     * @return The inverse matrix 
     */
    public final IMatrix inverseTriangleUpper(){
      if (this.numRow!=this.numCol){
          throw (new MismatchedSizeException("IMatrix.inverseTriangularSuperior() demands a square matrix."));
      }
      int n=this.numRow;
      IMatrix inverse=new IMatrix(n,n);
      for (int i=n-1;i>-1;i--){
        for (int j=n-1;j>-1;j--){
          if (i==j) inverse.setElement(i,j,1/this.matrix[i][i]);
          else if (i>j) inverse.setElement(i,j,0);
          else if (i<j) {
            double soma=0;
            for (int k=i+1;k<j+1;k++) soma+=inverse.getElement(k,j)*this.matrix[i][k];
            inverse.setElement(i,j,-soma/this.matrix[i][i]);
          }
        }
      }
      return inverse;
    }
    
    /**LU Decomposition for generic method with lii defined
     * This matrix must be a square matrix
     * The parameters matL and matU must be the same size of this matrix.     
     * @param matL The matrix into which the lower decomposition will be placed.
     * @param matU The matrix into which the upper decomposition will be placed.
     * @param c matL's value of the diagonal line 
     */
    public final void genericL(IMatrix matL, IMatrix matU, double c){
      if (this.numRow!=this.numCol){
          throw (new MismatchedSizeException("IMatrix.genericL() demands a square matrix."));
      }
      if (this.numRow!=matL.getNumRow() || this.numCol!=matL.getNumCol()) {
          throw (new MismatchedSizeException());
      }
      if (this.numRow!=matU.getNumRow() || this.numCol!=matU.getNumCol()) {
          throw (new MismatchedSizeException());
      }        
      int n=this.numRow;      
      for (int i=0;i<n;i++)matL.setElement(i,i,c);      
      double soma;
      for (int k=0;k<n;k++){ 
        for (int i=0;i<k;i++) {
          soma=0;        
          for (int j=0;j<i;j++) soma+=matL.getElement(k,j)*matU.getElement(j,i);
          matL.setElement(k,i,(this.matrix[k][i]-soma)/matU.getElement(i,i));
        }        
        for (int j=k;j<n;j++){
          soma=0;
          for (int i=0;i<k;i++) soma+=matL.getElement(k,i)*matU.getElement(i,j);
          matU.setElement(k,j,(this.matrix[k][j]-soma)/matL.getElement(j,j));          
        }
      }        
    }
    
    /**LU Decomposition for generic method with uii defined
     * This matrix must be a square matrix
     * The parameters matL and matU must be the same size of this matrix.     
     * @param matL The matrix into which the lower decomposition will be placed.
     * @param matU The matrix into which the upper decomposition will be placed.
     * @param c matU's value of the diagonal line 
     */
    public final void genericU(IMatrix matL, IMatrix matU, double c){
      if (this.numRow!=this.numCol){
          throw (new MismatchedSizeException("IMatrix.genericU() demands a square matrix."));
      }
      if (this.numRow!=matL.getNumRow() || this.numCol!=matL.getNumCol()) {
          throw (new MismatchedSizeException());
      }
      if (this.numRow!=matU.getNumRow() || this.numCol!=matU.getNumCol()) {
          throw (new MismatchedSizeException());
      }
      int n=this.numRow;      
      for (int i=0;i<n;i++)matU.setElement(i,i,c); 
      double soma;
      for (int k=0;k<n;k++){
        for (int i=0;i<k+1;i++) {
          soma=0;
          for (int j=0;j<i;j++) soma+=matL.getElement(k,j)*matU.getElement(j,i);
          matL.setElement(k,i,(this.matrix[k][i]-soma)/matU.getElement(i,i));
        }
        for (int j=k+1;j<n;j++){
          soma=0;
          for (int i=0;i<k;i++) soma+=matL.getElement(k,i)*matU.getElement(i,j);
          matU.setElement(k,j,(this.matrix[k][j]-soma)/matL.getElement(k,k));
        }
      }
    }
    
    /** Expand a vector(6x1) to a tensor(3x3). 
     *@param a The IVector(6x1) 
     *@return The IMatrix(3x3).
     */
     public IMatrix vecTensor3D(IVector a) {
         IMatrix Tensor = new IMatrix(3,3);
         Tensor.setElement(0,0,a.getElement(0));
         Tensor.setElement(1,1,a.getElement(1));
         Tensor.setElement(2,2,a.getElement(2));
         Tensor.setElement(0,1,a.getElement(3));
         Tensor.setElement(1,0,a.getElement(3));
         Tensor.setElement(0,2,a.getElement(4));
         Tensor.setElement(2,0,a.getElement(4));
         Tensor.setElement(1,2,a.getElement(5));
         Tensor.setElement(2,1,a.getElement(5));
         return Tensor;
     }
     
     public String toString(){
    	 DecimalFormat exp = new DecimalFormat ("0.00E00");
    	 String result = "";
    	 result += "\n";
    	 for (int i=0; i<numRow; i++){
    		 for (int j=0; j<numCol; j++){
    			 result += "  "+exp.format(matrix[i][j]);
    		 }
    		 result += "\n";	 
    	 }	 
    	 return result;
     }
	 
	 /** Returns the inverse of any square matrix
	     * @return The inverse matrix of this matrix
	     */
	 public final IMatrix pseudoInverse(){
	      IMatrix inv=new IMatrix(this.numCol,this.numRow);      
	      IMatrix AT=new IMatrix(this.numCol,this.numRow);
	      IMatrix ATA=new IMatrix(this.numCol,this.numCol);
		  AT.transpose(this);
		  ATA.mul(AT,this);
		  inv.mul(ATA.inverse(),AT);
	      return inv;
	 }
}
